      subroutine poisson(ncells, ijkcell, nvtxkm, nvtx, ijkvtx, itdim, iwid, &
         jwid, kwid, precon, tiny, relax, ibp1, jbp1, kbp1, cdlt, sdlt, strait&
         ,dx,dy, dz, c1x, c2x, c3x, c4x, c5x, c6x, c7x, c8x, c1y, c2y, c3y, c4y, c5y&
         , c6y, c7y, c8y, c1z, c2z, c3z, c4z, c5z, c6z, c7z, c8z, vol, vvol, &
         vvolaxis, q, qtilde, aqtilde, ht, itmax, error, srce, rdt, scalsq, wk&
         , phibc, ex, ey, ez, dive, residu, aq, phi, diag, wkl, wkr, wkf, wke, &
         periodicx)
!-----------------------------------------------
!   M o d u l e s
!-----------------------------------------------
      USE vast_kind_param, ONLY:  double
      use modify_com_M

!...Translated by Pacific-Sierra Research 77to90  4.3E  14:13:36   8/20/02
!...Switches: -yf -x1
      implicit none
!-----------------------------------------------
!   D u m m y   A r g u m e n t s
!-----------------------------------------------
      integer  :: ncells
      integer  :: nvtxkm
      integer  :: nvtx
      integer  :: itdim
      integer  :: iwid
      integer  :: jwid
      integer  :: kwid
      integer  :: ibp1
      integer  :: jbp1
      integer  :: kbp1
      integer , intent(in) :: itmax
      real(double)  :: tiny
      real(double) , intent(in) :: relax
      real(double)  :: cdlt
      real(double)  :: sdlt
      real(double)  :: strait
      real(double)  :: dx,dy,dz
      real(double) , intent(in) :: error
      real(double)  :: rdt
      real(double)  :: scalsq
      real(double)  :: wk
      real(double)  :: phibc
      real(double)  :: wkl
      real(double)  :: wkr
      real(double)  :: wkf
      real(double)  :: wke
      logical , intent(in) :: precon
      logical  :: periodicx
      integer  :: ijkcell(*)
      integer  :: ijkvtx(*)
      real(double)  :: c1x(*)
      real(double)  :: c2x(*)
      real(double)  :: c3x(*)
      real(double)  :: c4x(*)
      real(double)  :: c5x(*)
      real(double)  :: c6x(*)
      real(double)  :: c7x(*)
      real(double)  :: c8x(*)
      real(double)  :: c1y(*)
      real(double)  :: c2y(*)
      real(double)  :: c3y(*)
      real(double)  :: c4y(*)
      real(double)  :: c5y(*)
      real(double)  :: c6y(*)
      real(double)  :: c7y(*)
      real(double)  :: c8y(*)
      real(double)  :: c1z(*)
      real(double)  :: c2z(*)
      real(double)  :: c3z(*)
      real(double)  :: c4z(*)
      real(double)  :: c5z(*)
      real(double)  :: c6z(*)
      real(double)  :: c7z(*)
      real(double)  :: c8z(*)
      real(double)  :: vol(*)
      real(double)  :: vvol(*)
      real(double)  :: vvolaxis(*)
      real(double) , intent(inout) :: q(*)
      real(double)  :: qtilde(*)
      real(double)  :: aqtilde(*)
      real(double)  :: ht(*)
      real(double)  :: srce(*)
      real(double)  :: ex(*)
      real(double)  :: ey(*)
      real(double)  :: ez(*)
      real(double)  :: dive(*)
      real(double)  :: residu(*)
      real(double) , intent(inout) :: aq(*)
      real(double)  :: phi(*)
      real(double)  :: diag(*)
!-----------------------------------------------
!   L o c a l   V a r i a b l e s
!-----------------------------------------------
      integer :: nstart, n, k, j, i, ijk, iter, ijkmax
      real(double) :: lambda, rsum, bottom, rnorm, dxnorm, xnorm, alfa, resmax
!-----------------------------------------------
!
!
!
      nstart = (ibp1 - 1)*(jbp1 - 1) + 1
!
      if (precon) then
!
!
         call axisvol (ibp1, jbp1, iwid, jwid, kwid, vvol, nvtx, ijkvtx, vol, &
            vvolaxis)
!
         call diagonal (ncells, ijkcell, tiny, nstart, iwid, jwid, kwid, c1x, &
            c2x, c3x, c4x, c5x, c6x, c7x, c8x, c1y, c2y, c3y, c4y, c5y, c6y, &
            c7y, c8y, c1z, c2z, c3z, c4z, c5z, c6z, c7z, c8z, vvol, vol, diag)
!     &     c1z,c2z,c3z,c4z,c5z,c6z,c7z,c8z,vvolaxis,vol,diag)
!
      else
!
!     TEST CONJUGATE GRADIENT ITERATION
!
         diag(ijkcell(:ncells)) = 1.
!
!
      endif
!
!
      do k = 1, kbp1 + 1
         do j = 1, jbp1 + 1
!
            q(1+(j-1)*jwid+(k-1)*kwid:ibp1*iwid+1+(j-1)*jwid+(k-1)*kwid:iwid)&
                = 0.0
            aq(1+(j-1)*jwid+(k-1)*kwid:ibp1*iwid+1+(j-1)*jwid+(k-1)*kwid:iwid)&
                = 0.0
            qtilde(1+(j-1)*jwid+(k-1)*kwid:ibp1*iwid+1+(j-1)*jwid+(k-1)*kwid:&
               iwid) = 0.0
!      phi(ijk)=0.0
!
         end do
      end do
!
      call bcphi (ibp1, jbp1, kbp1, iwid, jwid, kwid, wkl, wkr, wkf, wke, phibc&
         , phi, periodicx)
!
!     CALCULATE THE INITIAL RESIDUAL ERROR
!
      call residue (ncells, ijkcell, nvtx, ijkvtx, iwid, jwid, kwid, ibp1, jbp1&
         , tiny, c1x, c2x, c3x, c4x, c5x, c6x, c7x, c8x, c1y, c2y, c3y, c4y, &
         c5y, c6y, c7y, c8y, c1z, c2z, c3z, c4z, c5z, c6z, c7z, c8z, vol, phi, &
         ex, ey, ez, vvol, dive, srce, rdt, scalsq, residu)
!
!
      rsum = 0.0
!
      do n = 1, ncells
         ijk = ijkcell(n)
         rsum = rsum + abs(residu(ijk))
      end do
!
      if (rsum < error) go to 2
!
!     ****************************************************************
!
!     begin conjugate gradient iteration
!
!     ****************************************************************
!
      do iter = 1, itmax
!
!     apply pre-conditioner
!
         do n = 1, ncells
!
!     sign is an experiment... normal condition _
!     is positive
!
            qtilde(ijkcell(n)) = residu(ijkcell(n))/diag(ijkcell(n))
!
            qtilde(ijkcell(n)) = phi(ijkcell(n)) + relax*qtilde(ijkcell(n))
!
         end do
!
         call bcphi (ibp1, jbp1, kbp1, iwid, jwid, kwid, wkl, wkr, wkf, wke, &
            phibc, qtilde, periodicx)
!
!
!
!     calculate residue
!
         call residue (ncells, ijkcell, nvtx, ijkvtx, iwid, jwid, kwid, ibp1, &
            jbp1, tiny, c1x, c2x, c3x, c4x, c5x, c6x, c7x, c8x, c1y, c2y, c3y, &
            c4y, c5y, c6y, c7y, c8y, c1z, c2z, c3z, c4z, c5z, c6z, c7z, c8z, &
            vol, qtilde, ex, ey, ez, vvol, dive, srce, rdt, scalsq, aqtilde)
!
!
         do n = 1, ncells
!
            aqtilde(ijkcell(n)) = residu(ijkcell(n)) - aqtilde(ijkcell(n))
!
            qtilde(ijkcell(n)) = qtilde(ijkcell(n)) - phi(ijkcell(n))
!
         end do
!
!     calculate LAMBDA
!
         bottom = 0.0
!
!
         lambda = sum(aqtilde(ijkcell(:ncells))*aq(ijkcell(:ncells))/diag(&
            ijkcell(:ncells)))
         bottom = sum(aq(ijkcell(:ncells))**2/diag(ijkcell(:ncells)))
!
!
         lambda = lambda/(bottom + 1.E-10)
!
!dir$ ivdep
         do n = 1, ncells
!
            q(ijkcell(n)) = qtilde(ijkcell(n)) - lambda*q(ijkcell(n))
            aq(ijkcell(n)) = aqtilde(ijkcell(n)) - lambda*aq(ijkcell(n))
!
         end do
!
!
         rnorm = 0.0
         dxnorm = 0.0
         xnorm = 0.0
         alfa = 0.0
         bottom = 0.0
!
!
         do n = 1, ncells
!
            ijk = ijkcell(n)
!
            rnorm = rnorm + abs(residu(ijkcell(n)))
            dxnorm = dxnorm + abs(qtilde(ijkcell(n))*lambda)
            xnorm = xnorm + abs(phi(ijkcell(n)))
            alfa = alfa + residu(ijk)*aq(ijk)/diag(ijk)
!
            bottom = bottom + aq(ijk)**2/diag(ijk)
!
         end do
!
         alfa = alfa/(bottom + 1.E-10)
!
!
!      write(*,*)iter,dxnorm,xnorm,error
!      write(*,*)iter,rnorm,error,rsum
         if (dxnorm<=error*xnorm .and. rnorm<=error*rsum) go to 2
!

         resmax = -1.E20
         ijkmax = 1
!
         do n = 1, ncells
            ijk = ijkcell(n)
            phi(ijk) = phi(ijk) + alfa*q(ijk)
            residu(ijk) = residu(ijk) - alfa*aq(ijk)
!
            if (resmax >= abs(residu(ijk))) cycle
            resmax = abs(residu(ijk))
            ijkmax = ijk
         end do
!
!
      end do
!     iteration failed
      write (*, *) 'poisson:  iteration relies on the good luck'
!
      do k = 1, kbp1 + 1
         do j = 1, jbp1 + 1
!
!            phi(1+(j-1)*jwid+(k-1)*kwid:ibp1*iwid+1+(j-1)*jwid+(k-1)*kwid:iwid)&
!                = 0.0
!
         end do
      end do
      return
!
    2 continue
      write (*, *) 'poisson:  iteration converged, iter=', iter
!
!     correct error in the divergence of e0
!
      call bcphi (ibp1, jbp1, kbp1, iwid, jwid, kwid, wkl, wkr, wkf, wke, phibc&
         , phi, periodicx)
!
!
!
      return
      end subroutine poisson
